library(SDMTools)
wd='/home/jc165798/working/Rob/';setwd(wd)

#read in the asciis
bdall=read.asc('Bdall.asc')
bdwet=read.asc('Bdwet.asc')
bc12=read.asc('bc12.asc')
bc01=read.asc('bc01.asc')
veg=read.asc('veg.asc')


######################################
#create the data
pos = as.data.frame(which(is.finite(bdall), arr.ind=TRUE))
pos$lat = getXYcoords(bdall)$y[pos$col]; pos$lon = getXYcoords(bdall)$x[pos$row] #append the lat lon
pos$bdall= extract.data(cbind(pos$lon,pos$lat),bdall) 	#extract the data
pos$bdwet= extract.data(cbind(pos$lon,pos$lat),bdwet)
pos$bc01= extract.data(cbind(pos$lon,pos$lat),bc01) 
pos$bc12 = extract.data(cbind(pos$lon,pos$lat),bc12)
pos$veg = extract.data(cbind(pos$lon,pos$lat),veg)


#read in the threshold to create sigdiff
threshold=read.csv('maxentResults.csv')
threshold=threshold$Minimum.training.presence.logistic.threshold
bdall[which(bdall<threshold[1])]=0; #set values below threshold to 0
bdwet[which(bdwet<threshold[2])]=0; #set values below threshold to 0
sig.diff=SigDiff(bdall,bdwet)
pos$sig.diff=extract.data(cbind(pos$lon,pos$lat),sig.diff)
write.csv(pos,'coord.csv')

######################################
#Create the plot
#edit and simplify the data
tdata = cbind(pos$bc01,pos$bc12); #
tdata[,1] = round(tdata[,1],1); tdata[,2] = round(tdata[,2]/25)*25 #round bc01 and bc12
tdata = unique(tdata) #find only unique combinations of rounded bc01 and bc12

pos$bdall[which(pos$bdall<threshold[1])]=0; #set values below threshold to 0
pos$bdwet[which(pos$bdwet<threshold[2])]=0; 
pos$bdall[which(pos$bdall>0)]=1 #set values above threshold to 1
pos$bdwet[which(pos$bdwet>0)]=1
tdata.all = unique(cbind(pos$bc01[which(pos$bdall==1)],pos$bc12[which(pos$bdall==1)]))#find only unique combinations
tdata.wet = unique(cbind(pos$bc01[which(pos$bdwet==1)],pos$bc12[which(pos$bdwet==1)]))#find only unique combinations


#make the plot
cols=colorRampPalette(c("tan","forestgreen"))(3)

png('compare.png')
	plot(tdata[,1],tdata[,2], xlab='Annual mean temperature', ylab='Annual rainfall', xlim=range(tdata[,1],na.rm=T),ylim=range(tdata[,2],na.rm=T), type='n', cex.lab=1, cex.axis=1, font.lab='2')
	points(tdata[,1],tdata[,2], col=cols[1])

	points(tdata.all,col=cols[2])

	points(tdata.wet,col=cols[3])
dev.off()

######################################
#make the image

out.range=range(sig.diff, na.rm=T)
ctan=out.range[2]*0.05
cgrey=out.range[2]*0.90
cgreen=out.range[2]*0.05

Colormap=c("grey",colorRampPalette(c("tan","forestgreen"))(100))
Colormap2=c(colorRampPalette(c("tan", 'burlywood'))(ctan), colorRampPalette(c("grey",'grey80'))(cgrey), colorRampPalette(c("forestgreen",'darkgreen'))(cgreen))	
pnts=cbind(x=c(429509,446345,446345,429509), y=c(8205059,8205059,8128594,8128594)) 

png ('sigdif_bd.png',width=3*dim(bdall)[1], height=1*dim(bdall)[2], bg= "white", pointsize=100)
    par(mfrow=c(1,3))
	  
	image(bdall, col=Colormap, ,xlab="", ylab="")
	legend.gradient(pnts,col=Colormap, c(0,1), title='Original')
	image(bdwet, col=Colormap, ,xlab="", ylab="")
	legend.gradient(pnts,col=Colormap, c(0,1), title='New')
	
	image(sig.diff,col=Colormap2, xlab="", ylab="")
	legend.gradient(pnts,col=Colormap2, c(1,0), title='SIG.DIFF')
dev.off()


